home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / CFIT4.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  2KB  |  106 lines
  1. program cfit4;        {164}
  2.  
  3. { plot service included }
  4. { Pascal program to perform a linear least-squares fit }
  5.  
  6. const    max    = 20;
  7.  
  8. type    index    = 1..max;
  9.     ary    = array[index] of real;
  10.  
  11. var    x,y,y_calc    : ary;
  12.     n        : integer;
  13.     first,done    : boolean;
  14.     a,b,correl_coef,
  15.     sigma_a,sigma_b,
  16.     see,seed    : real;
  17.  
  18. external procedure cls;
  19.  
  20. function random(dummy: integer): real;
  21. { random number 0-1 }
  22. { define seed=4.0 as global }
  23.  
  24. const    pi    = 3.14159;
  25.  
  26. var    x    : real;
  27.     i    : integer;
  28.  
  29. begin    { RANDOM }
  30.   x:=seed+pi;
  31.   x:=exp(5.0*ln(x));
  32.   seed:=x-trunc(x);
  33.   random:=seed
  34. end;    { RANDOM }
  35.  
  36.  
  37.  
  38. procedure get_data(var x,y: ary;
  39.            var n: integer);
  40. { get values for n and arrays x,y }
  41. { y is randomly scattered about a straight line }
  42.  
  43. const    a = 2.0;
  44.     b = 5.0;
  45.  
  46. var    i,j    : integer;
  47.     fudge    : real;
  48.  
  49. begin
  50.   write('Fudge? ');
  51.   readln(fudge);
  52.   if fudge<0.0 then done:=true
  53.   else
  54.     begin
  55.       repeat
  56.     write('How many points? ');
  57.     readln(n)
  58.       until (n>2) and (n<=max);
  59.       if first then first:=false else cls;
  60.       for i:=1 to n do
  61.     begin
  62.       j:=n+1-i;
  63.       x[i]:=j;
  64.       y[i]:=(a+b*j)*(1.0+(2.0*random(0)-1.0)*fudge)
  65.       end    { for-loop }
  66.     end        { if }
  67. end;        { procedure get_data }
  68.  
  69.  
  70. procedure write_data;
  71. { print out the answers }
  72. var    i,j    : integer;
  73.  
  74. begin
  75.   writeln;
  76.   writeln('  I      X       Y      YCALC');
  77.   for i:=1 to n do
  78.     writeln(i:3,x[i]:8:1,y[i]:9:2,y_calc[i]:9:2);
  79.   writeln;
  80.   writeln('Intercept is ',a:8:3,', sigma is ',sigma_a:8:3);
  81.   writeln('    Slope is ',b:8:2,', sigma is ',sigma_b:8:3);
  82.   writeln;
  83.   writeln('Correlation coeffizient is ',correl_coef:7:4);
  84.   for j:=1 to 40 do for i:=1 to 10000 do; cls
  85. end;        { write_data }
  86.  
  87. {$I C:LINFIT2.LIB }
  88.  
  89. {$I PLOT.LIB }
  90.  
  91. begin    { MAIN program }
  92.   seed:=4.0;
  93.   cls;
  94.   first:=true;
  95.   done:=false;
  96.   repeat
  97.     get_data(x,y,n);
  98.     if not done then
  99.       begin
  100.     linfit(x,y,y_calc,a,b,n);
  101.     write_data;
  102.     plot(x,y,y_calc,n);
  103.     end
  104.   until done
  105. end.
  106.